library(leaflet)  #パッケージの読み込み
popup = c("Robin", "Jakub", "Jannes")
leaflet() %>%
  addProviderTiles("NASAGIBS.ViirsEarthAtNight2012") %>%  #表示する地図の種類を選択
  addMarkers(lng = c(-3, 23, 11),  #経度設定
             lat = c(52, 53, 49),  #緯度設定
             popup = popup)
library(sf)  #ベクトルデータ
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(raster)  #ラスターデータ
## Loading required package: sp
library(spData)  #地図データ
library(spDataLarge)
world  #spDataに含まれているデータセット
## Simple feature collection with 177 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 11
##    iso_a2 name_long continent region_un subregion type  area_km2     pop lifeExp
##    <chr>  <chr>     <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>
##  1 FJ     Fiji      Oceania   Oceania   Melanesia Sove…   1.93e4  8.86e5    70.0
##  2 TZ     Tanzania  Africa    Africa    Eastern … Sove…   9.33e5  5.22e7    64.2
##  3 EH     Western … Africa    Africa    Northern… Inde…   9.63e4 NA         NA  
##  4 CA     Canada    North Am… Americas  Northern… Sove…   1.00e7  3.55e7    82.0
##  5 US     United S… North Am… Americas  Northern… Coun…   9.51e6  3.19e8    78.8
##  6 KZ     Kazakhst… Asia      Asia      Central … Sove…   2.73e6  1.73e7    71.6
##  7 UZ     Uzbekist… Asia      Asia      Central … Sove…   4.61e5  3.08e7    71.0
##  8 PG     Papua Ne… Oceania   Oceania   Melanesia Sove…   4.65e5  7.76e6    65.2
##  9 ID     Indonesia Asia      Asia      South-Ea… Sove…   1.82e6  2.55e8    68.9
## 10 AR     Argentina South Am… Americas  South Am… Sove…   2.78e6  4.30e7    76.3
## # … with 167 more rows, and 2 more variables: gdpPercap <dbl>,
## #   geom <MULTIPOLYGON [°]>
plot(world)
## Warning: plotting the first 9 out of 10 attributes; use max.plot = 10 to plot
## all

world["lifeExp"]  #lifeExpの配列を選択
## Simple feature collection with 177 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 2
##    lifeExp                                                                  geom
##      <dbl>                                                    <MULTIPOLYGON [°]>
##  1    70.0 (((180 -16.06713, 180 -16.55522, 179.3641 -16.80135, 178.7251 -17.01…
##  2    64.2 (((33.90371 -0.95, 34.07262 -1.05982, 37.69869 -3.09699, 37.7669 -3.…
##  3    NA   (((-8.66559 27.65643, -8.665124 27.58948, -8.6844 27.39574, -8.68729…
##  4    82.0 (((-122.84 49, -122.9742 49.00254, -124.9102 49.98456, -125.6246 50.…
##  5    78.8 (((-122.84 49, -120 49, -117.0312 49, -116.0482 49, -113 49, -110.05…
##  6    71.6 (((87.35997 49.21498, 86.59878 48.54918, 85.76823 48.45575, 85.72048…
##  7    71.0 (((55.96819 41.30864, 55.92892 44.99586, 58.50313 45.5868, 58.68999 …
##  8    65.2 (((141.0002 -2.600151, 142.7352 -3.289153, 144.584 -3.861418, 145.27…
##  9    68.9 (((141.0002 -2.600151, 141.0171 -5.859022, 141.0339 -9.117893, 140.1…
## 10    76.3 (((-68.63401 -52.63637, -68.25 -53.1, -67.75 -53.85, -66.45 -54.45, …
## # … with 167 more rows
world %>% dplyr::select(lifeExp) #上と同じ
## Simple feature collection with 177 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 2
##    lifeExp                                                                  geom
##      <dbl>                                                    <MULTIPOLYGON [°]>
##  1    70.0 (((180 -16.06713, 180 -16.55522, 179.3641 -16.80135, 178.7251 -17.01…
##  2    64.2 (((33.90371 -0.95, 34.07262 -1.05982, 37.69869 -3.09699, 37.7669 -3.…
##  3    NA   (((-8.66559 27.65643, -8.665124 27.58948, -8.6844 27.39574, -8.68729…
##  4    82.0 (((-122.84 49, -122.9742 49.00254, -124.9102 49.98456, -125.6246 50.…
##  5    78.8 (((-122.84 49, -120 49, -117.0312 49, -116.0482 49, -113 49, -110.05…
##  6    71.6 (((87.35997 49.21498, 86.59878 48.54918, 85.76823 48.45575, 85.72048…
##  7    71.0 (((55.96819 41.30864, 55.92892 44.99586, 58.50313 45.5868, 58.68999 …
##  8    65.2 (((141.0002 -2.600151, 142.7352 -3.289153, 144.584 -3.861418, 145.27…
##  9    68.9 (((141.0002 -2.600151, 141.0171 -5.859022, 141.0339 -9.117893, 140.1…
## 10    76.3 (((-68.63401 -52.63637, -68.25 -53.1, -67.75 -53.85, -66.45 -54.45, …
## # … with 167 more rows
world %>% dplyr::select(lifeExp) %>% st_drop_geometry()  #lifeExpの中身だけを表示
## # A tibble: 177 x 1
##    lifeExp
##  *   <dbl>
##  1    70.0
##  2    64.2
##  3    NA  
##  4    82.0
##  5    78.8
##  6    71.6
##  7    71.0
##  8    65.2
##  9    68.9
## 10    76.3
## # … with 167 more rows
st_drop_geometry(world["lifeExp"]) #上と同じ
## # A tibble: 177 x 1
##    lifeExp
##  *   <dbl>
##  1    70.0
##  2    64.2
##  3    NA  
##  4    82.0
##  5    78.8
##  6    71.6
##  7    71.0
##  8    65.2
##  9    68.9
## 10    76.3
## # … with 167 more rows
class(world)
## [1] "sf"         "tbl_df"     "tbl"        "data.frame"
plot(world["lifeExp"]) #geometryで表示

world_asia = world[world$continent == "Asia", ]  #Asiaだけを選択
asia = st_union(world_asia)  #Asiaの形を切り取り
## although coordinates are longitude/latitude, st_union assumes that they are planar
#> although coordinates are longitude/latitude, st_union assumes that they are planar

world_asia
## Simple feature collection with 47 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 26.04335 ymin: -10.35999 xmax: 145.5431 ymax: 55.38525
## Geodetic CRS:  WGS 84
## # A tibble: 47 x 11
##    iso_a2 name_long  continent region_un subregion type  area_km2    pop lifeExp
##    <chr>  <chr>      <chr>     <chr>     <chr>     <chr>    <dbl>  <dbl>   <dbl>
##  1 KZ     Kazakhstan Asia      Asia      Central … Sove… 2729811. 1.73e7    71.6
##  2 UZ     Uzbekistan Asia      Asia      Central … Sove…  461410. 3.08e7    71.0
##  3 ID     Indonesia  Asia      Asia      South-Ea… Sove… 1819251. 2.55e8    68.9
##  4 TL     Timor-Les… Asia      Asia      South-Ea… Sove…   14715. 1.21e6    68.3
##  5 IL     Israel     Asia      Asia      Western … Coun…   22991. 8.22e6    82.2
##  6 LB     Lebanon    Asia      Asia      Western … Sove…   10099. 5.60e6    79.2
##  7 PS     Palestine  Asia      Asia      Western … Disp…    5037. 4.29e6    73.1
##  8 JO     Jordan     Asia      Asia      Western … Sove…   89207. 8.81e6    74.0
##  9 AE     United Ar… Asia      Asia      Western … Sove…   79881. 9.07e6    76.9
## 10 QA     Qatar      Asia      Asia      Western … Sove…   11328. 2.37e6    77.9
## # … with 37 more rows, and 2 more variables: gdpPercap <dbl>,
## #   geom <MULTIPOLYGON [°]>
asia
## Geometry set for 1 feature 
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 26.04335 ymin: -10.35999 xmax: 145.5431 ymax: 55.38525
## Geodetic CRS:  WGS 84
## MULTIPOLYGON (((120.295 -10.25865, 118.9678 -9....
plot(world_asia)
## Warning: plotting the first 9 out of 10 attributes; use max.plot = 10 to plot
## all

plot(asia)

plot(world["pop"], reset = FALSE)  #worldのpopを地図にプロット
plot(asia, add = TRUE, col="red")  #Asiaだけを赤でプロット

multipoint_matrix = rbind(c(5, 2), c(1, 3), c(3, 4), c(3, 2))  #位置を行列にして表示
multi_point <- st_multipoint(multipoint_matrix)

multipoint_matrix
##      [,1] [,2]
## [1,]    5    2
## [2,]    1    3
## [3,]    3    4
## [4,]    3    2
class(multi_point)
## [1] "XY"         "MULTIPOINT" "sfg"
multi_point
## MULTIPOINT ((5 2), (1 3), (3 4), (3 2))
lnd_point = st_point(c(0.1, 51.5))                 #ロンドンの座標
class(lnd_point)
## [1] "XY"    "POINT" "sfg"
lnd_geom = st_sfc(lnd_point, crs = 4326)           #crs = coordinate reference system
lnd_attrib = data.frame(                           # data.frame object
  name = "London",
  temperature = 25,
  date = as.Date("2017-06-21")
  )
lnd_attrib
##     name temperature       date
## 1 London          25 2017-06-21
lnd_sf = st_sf(lnd_attrib, geometry = lnd_geom)    # sf object
lnd_sf
## Simple feature collection with 1 feature and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 0.1 ymin: 51.5 xmax: 0.1 ymax: 51.5
## Geodetic CRS:  WGS 84
##     name temperature       date         geometry
## 1 London          25 2017-06-21 POINT (0.1 51.5)
raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
new_raster = raster(raster_filepath)  #ラスターレイヤーを作成
new_raster
## class      : RasterLayer 
## dimensions : 457, 465, 212505  (nrow, ncol, ncell)
## resolution : 0.0008333333, 0.0008333333  (x, y)
## extent     : -113.2396, -112.8521, 37.13208, 37.51292  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : /usr/local/lib/R/site-library/spDataLarge/raster/srtm.tif 
## names      : srtm 
## values     : 1024, 2892  (min, max)
plot(new_raster)

multi_raster_file = system.file("raster/landsat.tif", package = "spDataLarge")
multi_raster_file
## [1] "/usr/local/lib/R/site-library/spDataLarge/raster/landsat.tif"
r_brick = brick(multi_raster_file)  #複数のラスターレイヤーを作成
r_brick
## class      : RasterBrick 
## dimensions : 1428, 1128, 1610784, 4  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 301905, 335745, 4111245, 4154085  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs 
## source     : /usr/local/lib/R/site-library/spDataLarge/raster/landsat.tif 
## names      : landsat.1, landsat.2, landsat.3, landsat.4 
## min values :      7550,      6404,      5678,      5252 
## max values :     19071,     22051,     25780,     31961
plot(r_brick)

ndvi = (r_brick[[4]] - r_brick[[3]])/(r_brick[[4]] + r_brick[[3]])  #正規化植生指数
plot(ndvi)

library(sf)
library(raster)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(stringr) # for working with strings (pattern matching)
library(tidyr)   # for unite() and separate()
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:raster':
## 
##     extract
sel_area = world$area_km2 < 10000  #10000km2未満のデータを抽出
summary(sel_area) # a logical vector
##    Mode   FALSE    TRUE 
## logical     170       7
#>    Mode   FALSE    TRUE 
#> logical     170       7
small_countries = world[sel_area, ]
small_countries
## Simple feature collection with 7 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -67.24243 ymin: -16.59785 xmax: 167.8449 ymax: 50.12805
## Geodetic CRS:  WGS 84
## # A tibble: 7 x 11
##   iso_a2 name_long  continent region_un subregion type  area_km2     pop lifeExp
##   <chr>  <chr>      <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>
## 1 PR     Puerto Ri… North Am… Americas  Caribbean Depe…    9225. 3534874    79.4
## 2 PS     Palestine  Asia      Asia      Western … Disp…    5037. 4294682    73.1
## 3 VU     Vanuatu    Oceania   Oceania   Melanesia Sove…    7490.  258850    71.7
## 4 LU     Luxembourg Europe    Europe    Western … Sove…    2417.  556319    82.2
## 5 <NA>   Northern … Asia      Asia      Western … Sove…    3786.      NA    NA  
## 6 CY     Cyprus     Asia      Asia      Western … Sove…    6207. 1152309    80.2
## 7 TT     Trinidad … North Am… Americas  Caribbean Sove…    7738. 1354493    70.4
## # … with 2 more variables: gdpPercap <dbl>, geom <MULTIPOLYGON [°]>
plot(small_countries)

plot(small_countries["pop"])  #popだけプロット

world %>% filter(area_km2 < 10000) %>% dplyr::select(pop) %>% plot()  #上までのプロットと同じ
##GDPperCapitalがtop10の国をピックアップして地図化してみましょう。
world %>%  top_n(n = 10, wt = gdpPercap) %>% 
  dplyr::select(gdpPercap) %>% 
  plot()

##GDPperCapitalがworst10の国をピックアップして地図化してみましょう。

world %>%  top_n(n = -10, wt = gdpPercap) %>%   #nを-にすると下から数えられる
  dplyr::select(gdpPercap) %>% 
  plot()

##日本はGDPperCapitalで上から何番目でしょう。
#1.gdpPercapをsort
gdp_sort <- world %>% arrange(desc(gdpPercap)) %>%   #arrangeで昇順にsort。descで降順
  dplyr::select(gdpPercap) %>% 
  st_drop_geometry()

#2.日本のgdpPercapを取り出したい
jp_gdp <- world %>% filter(name_long == "Japan") %>% 
  dplyr::select(gdpPercap) %>% 
  st_drop_geometry() %>% 
  as.numeric()

#3.gdpPercapが上からいくつか調べればいい
which(gdp_sort==jp_gdp)
## [1] 22

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)
##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.